home *** CD-ROM | disk | FTP | other *** search
/ TPUG - Toronto PET Users Group / TPUG Users Group CD / TPUG Users Group CD.iso / PET / S-Super PET / (s)t2.d64 / RESAMPLE.FTN < prev    next >
Text File  |  2009-01-18  |  17KB  |  491 lines
  1. *                           R E S A M P L E . F T N
  2. *                           ***********************
  3.  
  4. *     NOTE: REMOVE COMMENTS BEFORE RUNNING.  DATA ARRAYS NEED THE SPACE.
  5.  
  6. *     For purposes of time series analysis, the number  of  data  points
  7. *     must be an exact power of 2, at least when Fast Fourier Transform,
  8. *     and related filter algorithms are used.  However,  a  time  series
  9. *     will generally not consist of exactly the correct number  of  data
  10. *     points, particularly if the data has been obtained  through  tests
  11. *     or experimentation.  It is common to simply truncate the data file
  12. *     or pad zeros to the end in order to obtain the required number  of
  13. *     points.  Neither of these methods is entirely satisfactory.
  14.  
  15. *     The present program accepts a time series data file as  input  and
  16. *     creates from it  a  new  file  containing  any  number  of  points
  17. *     desired, including powers of  2.   When  an  increased  number  of
  18. *     points is requested, direct spline interpolation is used to obtain
  19. *     the new data points which will lie between the  old  data  points.
  20. *     When a smaller number of data points is requested, The time series
  21. *     is first passed through a low pass filter before  resampling  with
  22. *     the spline interpolation.  The filter cut-off frequency is set  at
  23. *     less than half the new sample  frequency  to  avoid  any  aliasing
  24. *     difficulties.
  25.  
  26. *     The program begins by requesting the name of the source data file.
  27. *     It is assumed that the source file will have been formatted by the
  28. *     program 'dataform.ftn' or 'datout'.  The four  header  records  of
  29. *     this file will be read and displayed to show the  number  of  data
  30. *     points, the time origin and the sample interval.  A  request  will
  31. *     then be made for  the  number  of  data  points  desired  and  the
  32. *     destination  file  name.   The  new  sample  interval  and  cutoff
  33. *     frequency will be displayed and the processing will proceed.   The
  34. *     new sample interval dt' is calculated from the old sample interval'
  35. *     dt using the relation:
  36.  
  37. *          dt' = dt * ( n - 1 ) / ( n' - 1 )
  38.  
  39. *     where n is the number of data points in the source file and n'  is'
  40. *     the number of points in the new file.  The  cutoff  frequency   fc
  41. *     for the original data file is calculated by assuming it  is  equal
  42. *     to half the sample frequency.  Thus:
  43.  
  44. *          fc =  (1/2) * 1 / dt
  45.  
  46. *     and this is expressed in cycles per second if dt  is  measured  in
  47. *     seconds.  if n' is less than n,  a  low  pass  filter  is  applied'
  48. *     before interpolation.  The new cutoff frequency will then be:
  49.  
  50. *          fc' = (1/2) * 1 / dt'
  51.  
  52. *     When this is expressed in  the  usual  normalized  form  with  the
  53. *     sample interval taken as unity, and  the  frequency  expressed  in
  54. *     radians per sample interval then:
  55.  
  56. *          omega' = fc' * pi * dt
  57.  
  58. *                 = (1/2) * pi * ( n' - 1 ) / ( n - 1)'
  59.  
  60. *     P R O G R A M   S T R U C T U R E
  61. *     *********************************
  62.  
  63. *     The program roughly conforms to the following steps:
  64.  
  65. *          Request name of source data file
  66. *          open Source data file
  67. *          read source data file header
  68. *          Display header information
  69. *          read source data file into data array
  70. *          close source file
  71. *          Request destination file information
  72. *          Check that n' .gt. 2  and n' .ne. n
  73.  
  74. *          if  n' .LT. n  THEN'
  75. *               Calculate new file cutoff frequency and filter length
  76. *               Display destination file header information
  77. *               Extend data array at both ends
  78. *               Apply lowpass filter to data array
  79. *          else
  80. *               Display destination file information
  81. *          endif
  82.  
  83. *          Save the filtered data in a temporary file
  84. *          execute subroutine to calculate spline sigma array
  85. *          Re-read the filtered data from the temporary file
  86.  
  87. *          open output file
  88. *          write output file header records
  89. *          Initialize data record counter
  90. *          while  record counter .lt. next to last record
  91. *               Calculate interpolated data one record at a time
  92. *               write output record
  93. *          endwhile
  94. *          Pad last record with zeros if necessary
  95. *          write last Record
  96. *          close output file
  97.  
  98. *          end
  99.  
  100.  
  101. *                        L O W P A S S    F I L T E R
  102. *                        ****************************
  103.  
  104. *     This subroutine replaces the  series   x   of  length   n    by  a
  105. *     smoothed version obtained with a  low-pass  digital  filter.   The
  106. *     parameter  m  is the maximum length of  x , that is, the dimension
  107. *     assigned to  x  in the calling  program.   The  filter  length  is
  108. *     (2 * ns + 1)  terms and the filter  weights  represent  the  least
  109. *     squares approximation to the cutoff filter with  cutoff  frequency
  110. *     w .  The cutoff frequency  w  is expressed in radians per  sample.
  111. *     The filter coefficients are proportional to:
  112.  
  113. *               sin( i * w )    sin( 2 * pi * i / ( 2 * ns + 1) )
  114. *               -----------  *  --------------------------------
  115. *                   i * w          2 * pi * i / ( 2 * ns + 1 )
  116.  
  117. *     where the first term  is  simply  the  Fourier  Transform  of  the
  118. *     rectangular pass band with cutoff frequency  w ,  and  the  second
  119. *     term is a convergence factor to  limit  the  Gibb's  oscillations.'
  120. *     The coefficients are  normalized  so  that  their  sum  is  unity,
  121. *     providing unity gain for a constant input array, that is, at  zero
  122. *     frequency.  The transfer function for  this  filter  will  have  a
  123. *     transition band of frequency width given by  4 * pi / (2 * ns + 1)
  124. *     which suggests that a  large  value  of   ns   would  reduce  this
  125. *     transition band width.
  126.  
  127. *     to avoid loss of data due to the finite length of the filter,  the
  128. *     input  array  is  extended  symmetrically  at  both  ends  by   ns
  129. *     positions.  This approach is not perfect since it  will  introduce
  130. *     some phase distortion in the filtered results at both  ends.   The
  131. *     extended array occupies  n + 2 * ns   positions  and  the  maximum
  132. *     size allowed for  x  must be at least as  large  as  the  extended
  133. *     array.  The filtered result is left in the  lower   n   positions.
  134. *     The size of  ns  is a balance  between  the  desire  for  a  small
  135. *     frequency transition band and a desire  to  limit  the  distortion
  136. *     length  at  the  array  ends  as  well  as   to   accelerate   the
  137. *     computations.  A cost may be defined consisting of the sum of  the
  138. *     proportion of the frequency transition band to the pass band,  and
  139. *     the proportion of the distortion length to the array  length.   An
  140. *     "optimum" value of  ns  may then be obtained  by  minimizing  this
  141. *     cost.  Thus:
  142.  
  143. *                         4 * pi / ( 2 * ns + 1 )           2 * ns
  144. *         cost   =   A *  -----------------------   +   B * ------
  145. *                                  w                          n
  146.  
  147. *     where  A  and  B  are weights which may  be  used  to  adjust  the
  148. *     importance of each term.  The minimization results in  the  filter
  149. *     half length of:
  150.  
  151. *          ns    =    ( SQRT( (A / B) * 4 * pi * n / w ) - 1 ) / 2
  152.  
  153. *     and the smallest integer part of this expression is taken.
  154.  
  155. *     The particular filter chosen  here  is  a  two  sided,  symmetric,
  156. *     non-recursive filter.  This  ensures  that  phase  distortion  and
  157. *     delays will not be introduced, except  at  the  end  points.   The
  158. *     filtering could be accomplished more rapidly  if  a  single  sided
  159. *     non-recursive filter were employed since FFT methods could then be
  160. *     used.  However, the time delays introduced this way could make  it
  161. *     difficult to relate the filtered data to real time.   Consequently
  162. *     the two sided filter has been chosen.
  163.  
  164.  
  165. *                  S P L I N E    I N T E R P O L A T I O N
  166. *                  ****************************************
  167.  
  168. *     The spline  interpolation  subroutine  is  based  on  the  program
  169. *     published  in  the  book  "Computer   Methods   for   Mathematical"
  170. *     Computations" by Forsythe, Malcolm and Moler, Prentice-Hall  1977,"
  171. *     Page  76.   The  program  has  been  modified  to  reduce  storage
  172. *     requirements  by  capitalizing  on  the  fact  that  the  data  is
  173. *     available at equally spaced intervals.  Further more, rather  than
  174. *     calculating and storing the three cubic coefficients for each data
  175. *     point, only the so called 'sigma' array is stored.  That  is,  the
  176. *     interpolated points are obtained from the following  interpolation
  177. *     equation:
  178.  
  179. *         x(z)  =  w * x(i+1)  +  w' * x(i)'
  180.  
  181. *                    +  h*h * [w*(w*w-1)*s(i+1) + w'*(w'*w'-1)*s(i)]'
  182.  
  183. *     where
  184. *                 h  = sample interval
  185.  
  186. *                 w  = z - h * i
  187.  
  188. *                 w' = 1 - w'
  189.  
  190. *               s(i) = i'th sigma quantity'
  191.  
  192. * ************************************************************************
  193.  
  194.       real x(1700), s(1700), tzero, dt, t(5)
  195.       character source$file, destination$file, file$id, var$fmt, garbag
  196.  
  197. * DIMENSIONS OF  x  AND  s  AND THE PARAMETER  m   DETERMINE MAX ARRAY SIZE.
  198.  
  199.       m = 1700
  200.       pi = 4.0 * atan(1.0)
  201.  
  202. * CLEAR SCREEN, IDENTIFY PROGRAM AND REQUEST DATA FILE NAME.
  203.  
  204.       write(6,"'1'")
  205.       print,"        P R O G R A M    R E S A M P L E . F T N"
  206.       print,"        ****************************************"
  207.       print
  208.       print," Enter the source file name"
  209.       read, source$file
  210.  
  211. * READ THE INPUT DATA FILE.
  212.  
  213.       call data$input(source$file, file$id, m, var$fmt, tzero, dt, fc, x, n)
  214.  
  215. * HOW MANY POINTS IN RE-SAMPLED FILE?  ENSURE NUMBER IS IN ACCEPTABLE RANGE.
  216.  
  217.       loop
  218.           write(6,1) char(11),char(11)
  219.     1     format(1x,2a1,"Enter number of data points for the destination file")
  220.           read, nprime
  221.           quitif (nprime .gt. 2) .and. (nprime .ne. n)
  222.           write(6,"'+          '")
  223.       endloop
  224.  
  225. * ENTER AND DISPLAY OUTPUT FILE HEADER INFORMATION.
  226. * APPLY LOW PASS FILTER IF RESAMPLING AT A REDUCED RATE.
  227.  
  228.       print,"Enter destination file name"
  229.       read, destination$file
  230.       file$id = file$id//"- RESAMPLED"
  231.       dtprime = dt * float(n - 1) / float(nprime - 1)
  232.       print
  233.       print,"DESTINATION FILE PARAMETERS"
  234.       if nprime .lt. n
  235.           fcprime = 3600.0 * 0.5 / dtprime
  236.           call display(file$id, nprime, var$fmt, tzero, dtprime, fcprime)
  237.           w = 0.5 * pi * float(nprime-1) / float(n-1)
  238.           call lowpass(x, n, m, w)
  239.       else
  240.           fcprime = fc
  241.           call display(file$id, nprime, var$fmt, tzero, dtprime, fcprime)
  242.       endif
  243.  
  244. * CALCULATE SPLINE INTERPOLATION ARRAY FOR EQUALLY SPACED DATA BUT FIRST
  245. * SAVE THE FILTERED  x  ARRAY IN A TEMPORARY FILE SINCE splineq WILL
  246. * DESTROY IT.  RE-READ THE DATA WHEN splineq IS FINISHED.
  247.  
  248.       open(unit=8,file="temporary.dat")
  249.       write(8,var$fmt) (x(i),i=1,n)
  250.       close(unit=8)
  251.       call splineq(n, m, x, s)
  252.       open(unit=8,file="temporary.dat")
  253.       read(8,var$fmt) (x(i),i=1,n)
  254.       close(unit=8)
  255.  
  256. * WRITE THE FOUR OUTPUT FILE HEADER RECORDS.  USE THE STANDARD FORMATS.
  257.  
  258.       open(unit=8,file=destination$file)
  259.       write(8,"a79,'1'") file$id
  260.       write(8,"i5,74x,'2'") nprime
  261.       write(8,"a79,'3'") var$fmt
  262.       write(8,"3f10.5,49x,'4'") tzero, dtprime, fcprime
  263.  
  264. * OBTAIN INTERPOLATED POINTS AND WRITE TO FILE, RECORD BY RECORD.
  265.  
  266.       l = 4
  267.       factor = float(n-1)/float(nprime-1)
  268.       jend = nprime/5
  269.       do j = 1 , jend
  270.           k = 5 * (j - 1) - 1
  271.           do i = 1 , 5
  272.               u = 1.0 + float(k+i) * factor
  273.               t(i) = val(n, u, x, s)
  274.           enddo
  275.           l = j + 4
  276.           write(8,"5e15.8,i5") (t(i),i=1,5), l
  277.       enddo
  278.  
  279. * PAD LAST RECORD WITH ZEROS IF NEEDED.
  280.  
  281.       k = 5 * jend
  282.       if k .lt. nprime
  283.           k = k - 1
  284.           do i = 1 , 5
  285.               t(i) = 0.0
  286.               u = 1.0 + float(k+i) * factor
  287.               if (k+i) .lt. nprime
  288.                   t(i) = val(n, u, x, s)
  289.               endif
  290.           enddo
  291.           l = l + 1
  292.           write(8,"5e15.8,i5") (t(i),i=1,5), l
  293.       endif
  294.  
  295.       close(unit=8)
  296.       end
  297.  
  298.       subroutine data$input(file$name, file$id, m, var$fmt, tzero, dt, fc, x, n)
  299.       character file$name, file$id, var$fmt, two$characters, two$blanks
  300.       integer n, character$count
  301.       real tzero, dt, fc, x(m)
  302.       open(unit=8,file=file$name)
  303.  
  304. * READ FILE ID AND REMOVE TRAILING BLANKS.
  305.  
  306.       read(8,"a72") file$id
  307.       character$count = 0
  308.       two$blanks = "  "
  309.       while character$count .lt. 72
  310.           character$count = character$count + 1
  311.           two$characters = substr(file$id, character$count, 2)
  312.           quitif two$characters = two$blanks
  313.       endwhile
  314.       file$id = substr(file$id, 1, character$count)
  315.  
  316. * READ REMAINING THREE HEADER RECORDS.
  317.  
  318.       read(8,"i5") n
  319.       read(8,"a9") var$fmt
  320.       read(8,"3f10.5") tzero, dt, fc
  321.  
  322. * CALCULATE ASSUMED FILE CUTOFF FREQUENCY IF NOT RECORDED.
  323.  
  324.       if  fc = 0.0
  325.           fc = 3600.0 * 0.5 / dt
  326.       endif
  327.  
  328. * DISPLAY THE HEADER INFORMATION.
  329.  
  330.       call display(file$id, n, var$fmt, tzero, dt, fc)
  331.  
  332. * READ THE DATA.
  333.  
  334.       read(8,var$fmt) (x(j),j=1,n)
  335.       close(unit=8)
  336.       return
  337.       end
  338.  
  339.  
  340.       subroutine display(id,n,vform,t,dt,f)
  341.       character id,vform
  342.       print
  343.       print,"      File Identifier: ", id
  344.       print,"Number of Data Points: ", n
  345.       print,"          Data Format: ", vform
  346.       print,"          Time Origin: ", t
  347.       print,"            Time Step: ", dt
  348.       write(6,1) f
  349.     1 format("      Cutoff frequency: ", f6.2, " Cycles per Hour"///)
  350.       return
  351.       end
  352.  
  353.  
  354.       subroutine lowpass(x, n, m, w)
  355.       real x(m), h(100)
  356.       integer hsize
  357.       hsize = 100
  358.       pi = 4.0 * atan(1.0)
  359.  
  360. * CALCULATE FILTER HALF LENGTH.
  361.  
  362.       A = 1.0
  363.       B = 10.0
  364.       ns = ( sqrt( (A / B) * 4 * pi * float(n) / w ) -1 ) / 2.0
  365.       length = 2 * ns + 1
  366.       if  length .gt. hsize
  367.           print,"Filter length =",length," exceeds dimension specification"
  368.           stop
  369.       endif
  370.       if  (n + 2 * ns) .gt. m
  371.           print,"Extended  x  array exceeds available dimension"
  372.           stop
  373.       endif
  374.  
  375. *  EXTEND THE ARRAY.
  376.  
  377.       n1 = n+ns+1
  378.       n2 = n+1
  379.       do i = 1 , n
  380.           x(n1-i) = x(n2-i)
  381.       enddo
  382.       n1 = ns+1
  383.       n2 = n+ns
  384.       do i = 1 , ns
  385.           x(n1-i) = x(n1+i)
  386.           x(n2+i) = x(n2-i)
  387.       enddo
  388.  
  389. * CALCULATE FILTER WEIGHTS.
  390.  
  391.       hzero = w / pi
  392.       con = 2.0 * pi / float(2*ns+1)
  393.       sum = hzero
  394.       do i = 1 , ns
  395.           h(i) = hzero * sinc( float(i) * w ) * sinc( float(i) * con )
  396.           sum  = sum + 2.0 * h(i)
  397.       enddo
  398.       hzero = hzero / sum
  399.       do i = 1 , ns
  400.           h(i) = h(i) / sum
  401.       enddo
  402.  
  403. * APPLY FILTER WEIGHTS.
  404.  
  405.       do i = 1 , n
  406.           ins  = i + ns
  407.           temp = hzero * x(ins)
  408.           do j = 1 , ns
  409.               temp = temp + (x(ins+j) + x(ins-j)) * h(j)
  410.           enddo
  411.           x(i) = temp
  412.       enddo
  413.       return
  414.       end
  415.  
  416.       function sinc(z)
  417.           sinc = sin(z) / z
  418.       end
  419.  
  420.  
  421.  
  422.       subroutine splineq(n, m, x, s)
  423.       real x(m), s(m)
  424.  
  425. * NOTE: x  ARRAY IS DESTROYED BY THIS ROUTINE.
  426.  
  427.       if (n .le. 3)
  428.           print," Array length is less than or equal to 3.  This"
  429.           print," spline routine needs at least 4 points."
  430.           print," Execution has been terminated."
  431.           stop
  432.       endif
  433.  
  434. *  SETUP TRIDIAGONAL SYSTEM FOR EQUALLY SPACED INTERVALS.
  435.  
  436.       nm1 = n - 1
  437.       s(2) = x(2)-x(1)
  438.       bi = 4.0
  439.       do i = 2 , nm1
  440.           s(i+1) = x(i+1) - x(i)
  441.           s(i) = s(i+1) - s(i)
  442.       enddo
  443.  
  444. *  END CONDITIONS.
  445.  
  446.       b1 = -1.0
  447.       bn = -1.0
  448.       s(1) = ( s(3)   - s(2)   ) / 2.0
  449.       s(n) = ( s(n-1) - s(n-2) ) / 2.0
  450.       s(1) =  s(1) / 3.0
  451.       s(n) = -s(n) / 3.0
  452.  
  453. *  FORWARD ELIMINATION.  THIS IS WHERE THE  x  ARRAY IS USED FOR STORAGE
  454. *  OF THE DIAGONAL  b  ELEMENTS.
  455.  
  456.       b = b1
  457.       x(1) = b
  458.       do i = 2 , n
  459.           t = 1.0/b
  460.           b = bi - t
  461.           x(i) = b
  462.           s(i) = s(i)-t*s(i-1)
  463.       enddo
  464.       bn = bn - t
  465.  
  466. *  BACK SUBSTITUTION.
  467.  
  468.       s(n) = s(n)/bn
  469.       b = 1.0/t
  470.       do j = 1 , nm1
  471.           i = n-j
  472.           b = x(i)
  473.           s(i) = (s(i)-s(i+1))/b
  474.       enddo
  475.       return
  476.       end
  477.  
  478.       real function val(n, u, x, s)
  479.       real x(n), s(n), u
  480.       integer i
  481.       i = u
  482.       if (i .ge. n) i = n - 1
  483.       w = u - float(i)
  484.       wbar = 1.0 - w
  485.       temp = w * x(i+1) + wbar * x(i)
  486.       temp = temp + w*(w+1.0)*(w-1.0)*s(i+1)
  487.       temp = temp + wbar*(wbar+1.0)*(wbar-1.0)*s(i)
  488.       val  = temp
  489.       return
  490.       end
  491.